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Abstract: We develop a detailed and comprehensive description of neutrino oscillations 
^ . driven by the 1-3 mixing in the matter of the Earth. The description is valid for the realistic 

(PREM) Earth density profile in the whole range of nadir angles and for neutrino energies 
above 1 GeV. It can be applied to oscillations of atmospheric and accelerator neutrinos. The 
results are presented in the form of neutrino oscillograms of the Earth, i.e. the contours of 
equal oscillation probabilities in the neutrino energy-nadir angle plane. A detailed physics 
interpretation of the oscillograms, which includes the MSW peaks, parametric ridges, local 
maxima, zeros and saddle points, is given in terms of the amplitude and phase conditions. 
Precise analytic formulas for the probabilities are obtained. We study the dependence 
of the oscillation pattern on #13 and find, in particular, that the transition probability 
P > 1/2 appears for sin 2 2^13 as small as ~ 0.009. We consider the dependence of the 
oscillation pattern on the matter density profile and comment on the possibility of the 
oscillation tomography of the Earth. 
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1. Introduction 

Substantial future progress in neutrino physics will be related to the long baseline exper- 
iments as well as studies of the cosmic and atmospheric neutrinos. The key element of 
these experiments is that neutrinos propagate long distances inside the Earth before reach- 
ing detectors. Oscillations in the matter of the Earth change flavor properties of neutrino 
fluxes, which opens up possibilities to 
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• study dynamics of various oscillation effects; 

• determine yet unknown oscillation parameters: 1-3 mixing, deviation of the 2-3 mixing 
from maximal and its octant, the type of the neutrino mass hierarchy, and CP-violation; 

• perform the oscillation tomography of the Earth; 

• search for various new effects caused by non-standard neutrino interactions and/or by 
exotic physics (such as violation of Lorentz or CPT invariance) . 

The sensitivity of present-day experiments to the corresponding effects, that appear 
usually as small corrections to the leading oscillation phenomena (e.g. <-> v T vacuum 
oscillations), is rather low. In the case of atmospheric neutrinos this is related to low 
fluxes and therefore low statistics of events, and also to uncertainties in the fluxes. The 
accelerator neutrino beams are, in principle, well controlled, however in majority of the 
projects the baselines are relatively short, whereas the very long baseline experiments (e.g. 
with neutrino factories) employ high energy neutrinos. This means that the region of 
energy E v ~ (3 — 10) GeV and nadir angle Q v = — 70°, where the most interesting 
oscillation phenomena occur, is not covered. Therefore one needs to rely on high precision 
measurements of the suppressed "tails" of these phenomena. The dilemma is whether 
to study small effects with a number of systematical errors and degeneracies, or develop 
experimental approaches that will allow to cover the regions of large oscillation effects. The 
latter option is realized in the concepts of very long baseline accelerator experiments and 
very large volume atmospheric neutrino detectors. 

In connection with possible future experimental developments and discussions of new 
strategies of research, a detailed and comprehensive study of physics of the neutrino oscil- 
lations in the Earth is needed. Some studies in this direction have been carried out in the 
past. 

In the matter of the Earth, the resonance enhancement of v e <-> v„ and v T oscillations 
can take place [1, 2], with the MSW resonance peak at E u /Am 2 ~ 2.5 xlO 3 GeV/eV 2 [3, 4]. 
It was also recognized that the size of the Earth is comparable with the neutrino refraction 
length, and consequently strong enhancement of oscillations can occur for rather deep 
trajectories (small nadir angles) and not too small mixing angles: sin 2 26^ > 0.03. With 
decreasing baseline, first the matter effect and then the oscillation effect disappear. This 
phenomenon [1, 5], termed "vacuum mimicking", implies that for short baselines the 
flavour transitions are described by the vacuum oscillation formula. 

For 1-2 mixing the resonance is in the neutrino channel and at rather low energies: E v ~ 
0.1 GeV. For non-zero 1-3 mixing the MSW resonance peak related to the atmospheric mass 
splitting is at E v ~ 6 GeV. The peak is narrow and, depending on the mass hierarchy, the 
enhancement occurs in the neutrino (normal hierarchy) or antineutrino (inverted hierarchy) 
channels. Possibilities to observe this resonance peak have been explored in connection 
to the measurements of the mixing angle #13 and determination of the neutrino mass 
hierarchy [6, 7, 8, 9, 10, 11, 12]. Furthermore, it was realized that matter effects can also 
strongly influence <-> v T oscillations [13, 14]. 
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A qualitatively new oscillation phenomenon can be realized for neutrinos crossing the 
core of the Earth due to a sharp change of the density at the border between the core and 
the mantle [15, 16, 17, 18, 19, 20]. In particular, for non-zero 1-3 mixing, the existence of an 
additional peak is predicted in the energy distribution between the peaks due to the MSW 
resonances in the core and mantle. This core-mantle effect was interpreted as being due 
to the parametric enhancement of neutrino oscillations [15, 17, 18]. 1 The additional peak 
is a manifestation of the parametric resonance for 3 layers (1.5 period of the "castle wall" 
profile) in matter [21]. The parametric resonance occurs when the variation of the matter 
density along the neutrino trajectory is in a certain way correlated with the values of the 
oscillation parameters [22, 23, 24]. A different interpretation of this peak, as being due to 
an interference between the contributions to the oscillation amplitude from different layers 
of the Earth's mantle and core, has been discussed in [16, 19, 20]. It was uncovered [25] 
that the parametric resonance condition is fulfilled also at large energies (above the mantle 
MSW resonance), leading to the appearance of two parametric peaks in the nadir angle 
distribution. For general consideration of evolution in the multi-layer media, see [26]. 

Apart from the large scale structures of the density profile, effects of small scale density 
perturbation have been explored [27]. 

Analytic approaches have been developed to describe the physics of various oscillation 
effects, to understand the influence of the density profile modifications on the oscillation 
probabilities and also to simplify numerical computations. Many studies have been per- 
formed in the constant density approximation [3, 4] or approximation of several layers of 
constant densities [15, 16, 17, 18, 19, 20, 23, 24, 28, 13, 29]. For the varying densities within 
the layers a perturbation-theory approach has been developed in [30] for the description of 
the solar neutrino oscillations in the Earth. 

For 3is mixing the analytic approaches employ various expansions in the small param- 
eters sin 2 6>i3 and/or in the ratio of the mass squared differences Am^/Am^ [31]. For the 
case of non-constant density matter, different perturbation theory approaches have been 
developed for oscillations in low density [32, 33] and high density [34, 25] media. In [25] 
an analytic description of oscillations in the high energy limit has been worked out for the 
case of the realistic (PREM) density profile. The influence of the Earth density profile on 
the oscillation probabilities has also been considered in [35]. 

To some extent, the results obtained so far had fragmentary character. This paper 
is the first one in a series of papers we devote to a detailed and comprehensive study 
of oscillations of neutrinos inside the Earth. Here we present a thorough description of 
neutrino oscillations caused by non-zero 1-3 mixing, with the emphasis on the physics of 
the phenomenon. In particular, we study the complex pattern and interplay of various 
oscillation resonances. We show that the main features of this complex pattern can be 
easily understood in terms of different realizations of just two simple conditions - the am- 
plitude condition and the phase condition. We perform our study in terms of "neutrino 
oscillograms" of the Earth: contours of constant oscillation probabilities in the plane of 

lr The parametric resonance in neutrino oscillations in the Earth was first discussed for <-> v s oscil- 
lations of atmospheric neutrinos in [15], though the parametric peak appeared in some early numerical 
computations [3, 4]. 
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neutrino energy and nadir angle. These plots have been presented in several earlier pub- 
lications as illustrations [19, 8, 36]. Here we use them as the main tool of the study. We 
present a description of all the structures of the oscillograms and of their dependence on 
013 and the density profile. We generalize the amplitude and the phase conditions to the 
case of non-constant densities in the layers of matter and develop an appropriate analytic 
formalism. 

The paper is organized as follows. In Sec. 2, after some generalities, we introduce 
the neutrino oscillograms of the Earth, concentrating mainly on u e <-> v ii {y T ) transitions. 
We discuss the accuracy of the constant-density layers approximation of the Earth matter 
profile and describe the graphical representation of the conversion. In Sec. 3 we give 
the physics interpretation of the oscillograms. In Sec. 4 we derive approximate analytic 
formulas for probabilities for a realistic matter density profile. In Sec. 5 we study the 
dependence of the oscillograms on the 1-3 mixing and the density profile of the Earth, and 
we discuss the oscillation probabilities for the other channels. Discussion and conclusions 
follow in Sec. 6. 



2. Neutrino oscillograms of the Earth 
2.1 Context. Evolution matrix 

We consider the three-flavor neutrino system with the state vector Vj = (u e , v T ) T . Its 
evolution in matter is described by the equation 

.dvt /UM 2 W A , n ,. 

'-d = {-^r +y ^r" <"> 

where E v is the neutrino energy and M 2 = diag(0, Am^, Am^) is the diagonal matrix of 
neutrino mass squared differences. V(x) = diag(V(x), 0, 0) is the matrix of matter-induced 
neutrino potentials with V(x) = y/2GpN e (x), Gf and N e (x) being the Fermi constant and 
the electron number density, respectively. The mixing matrix U defined through Vf = Uu m , 
where v m = (yi, 1/2, v?,) T is the vector of neutrino mass eigenstates, can be parameterized 
as 

U = U 2 3lsU 13 I^U 12 . (2.2) 

Here the matrices Uij = Uij(9ij) describe rotations in the ij-planes by the angles 6ij, and 
Is = diag(l, 1, e tS ), where 5 is the Dirac-type CP-violating phase. 

Let us introduce the evolution matrix (the matrix of transition and survival amplitudes) 
S(x,xo), which describes the evolution of the neutrino state over a finite distance: from 
xq to x. To simplify the presentation, throughout the paper we will use the notation 
S(x) = S(x,0) and S = S(L), where L is the total length of the trajectory. The matrix 
S(x) satisfies the same evolution equation as the state vector (2.1): 

i^M. = H(x) S(x) . (2.3) 
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The solution of equation (2.3) with the initial condition 5(0) = 1 can be formally written 



as 

(•X 



5(x) =Texp l-i J Hdxj . (2.4) 

It is convenient to consider the evolution of the neutrino system in the propagation 
basis, v = (is e , 9 2 , v 3 ) r ', defined through the relation 

vf = U 23 I s v • (2.5) 

As follows from (2.1) and (2.2), the Hamiltonian H, that describes the evolution of the 
neutrino vector of state z>, is 

H(x) = ^-U 13 U 12 M 2 ul 2 Ul 3 + V(x) . (2.6) 

This Hamiltonian does not depend on the 2-3 mixing and CP-violating phase. The depen- 
dence on these parameters appears when one projects the initial state on the propagation 
basis and the final state back onto the original flavor basis. Explicitly, the Hamiltonian H 
reads 

/ s i3 + s i2 c i3 r A + 2V(x)E„/Am 2 , 1 su c 12 c 13 r A s 13 ci 3 (l - s\ 2 r A )\ 



2E U 



c\ 2 r A -S12 Cl2 si 3 r A 

V • • • • • • C 13 + S 12 s 13 r A J 



(2.7) 



where Cij = cos 6ij , Sij = sin 6-ij and 



_ Am| 

rA = • (2 ' 8) 

According to (2.5), the evolution matrix S(x) in the basis (y e ,V2,vz) is related to S(x) by 
the transformation: 

S(x) = US(x)U\ U = U 23 Is- (2.9) 

The evolution of S(x) is given by the equation analogous to Eq. (2.3) with the Hamiltonian 
H(x). 

2.2 High energy neutrino approximation 

For sufficiently high energies (E v > 1 — 2 GeV) one can neglect the 1-2 mass splitting. 
Then, according to (2.7), the state v 2 decouples from the rest of the system and does not 
evolve. Therefore, if we parameterize 5 in the basis (is e , 9 2 , v 3 ) as 

' A ee A e2 A e3 ^ 

A 2e A 22 A 23 (2.10) 
K A 3e A 32 A 33 ) 

we find that in this approximation A e2 = A 2e = A 23 = A 32 = 0, A 22 = 1 and the evolution 
matrix in the flavor basis takes the form 

( A ee s 23 A e3 c 23 A e3 

5» S 23 ^3e c 2 23 A 22 + s 2 23 A 33 -s 23 c 23 (A 22 - A 33 ) I , (2.11) 
\c23-43e s 23 c 23 (A 22 - A 33 ) s 2 23 A 22 + 4 3 A 33 
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where we omitted the CP-phase factor e lS since CP-violating effects are absent in the limit 
r a — > 0. We will consider the complete 3v system in the next publication [37]. Denoting 



P A = \A e3 \ 2 = \A 3e \ 2 , (2.12) 



we obtain from (2.11) 





-+ v e ) 


= 1 - Pa, 


P{Ve~ 


+ v v) 


= P(^ ^e) = S23P4 , 


P{Ve~ 




= P(v T -»■ f e ) = C23P4 , 


P(^- 


+ v v) 


= l-4P4-244 [l-ReA 33 ] 




-> I/ r ) 


= -^L 4$ P4 + 2s 23C23 [1 - Re A 33 



(2.13) 
(2.14) 
(2.15) 
(2.16) 
(2.17) 

The formulas in Eqs. (2.13-2.17) reproduce the probabilities derived in [18]. 

Thus, in the approximation Am^ = the dynamical problem is reduced to two 
flavor evolution problem. Throughout this paper we work in the basis where, for a 2-flavor 
neutrino system, H\\ = — P22 = — cos 29isAm 2 /4E + V/2. This can be always achieved by 
subtracting from H a matrix proportional to the unit matrix and correspondingly rephasing 
the neutrino vector of state. This symmetric Hamiltonian is related to the Hamiltonian of 
the 1-3 subsystem H( 13 \ obtained from Eq. (2.7) by taking the limit = and removing 
the decoupled state z>2, as 

For the 2v case the unitary evolution matrix can be parameterized as 

\a\ 2 + \(3\ 2 = 1. (2.19) 




For density profiles that are symmetric with respect to the midpoint of the neutrino trajec- 
tory (for brevity, symmetric profiles) , T invariance leads to the equality of the off-diagonal 
elements of the evolution matrix, (3 = —(3*, which means that (3 is pure imaginary [38]. 
For a single layer of constant density the solution can be written explicitly: 

_,. . / cos6(x) + icos29 m sm6(x) — ism29 m smd(x) \ , , . 
b (x) = , <p(x) = to x . 

V — ism29 m sm<f>(x) cos <fi(x) — icos29 m sm(f>{x) I 

(2.20) 

Here is 9 m is the mixing angle in matter and <f>(x) is the half-phase of oscillations in matter 
with 

° = ^ - fi^W~W^¥^Wf- <2 - 2i) 

The moduli squared of the elements of S reproduce the well-known probabilities for oscil- 
lations in a uniform medium. Thus in the notation of Eq. (2.19), 

a = cos 4> + i cos 29 m sin <fi , (3 = — i sin 29 m sin <f> . (2.22) 
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In what follows we will generalize this result to the cases of several layers of constant 
densities and also of changing densities within the layers. 

The scale of the matter effects in neutrino oscillations is set up by the refraction length, 
lo = 2it/(\/2Gf N e ), which is comparable with the radius of the Earth R = 6371 km. The 
oscillation length in matter l m is given by 

^ = 2w . (2.23) 

At high energies the length l m essentially coincides with the refraction length: l m ~ Iq. 
In the resonance channel, with decreasing energy l m increases and reaches its maximum 
slightly above the resonance energy. At the MSW resonance, E u = Er, one has l m = 
l v / sin 2#i3, where l v is the vacuum oscillation length. Below the resonance, with energy 
further decreasing, l m decreases and approaches l v : l m (E v <C Er) ps l v . 

In a non-uniform density medium the adiabatic half-phase of oscillations is defined as 

<p{x)= uj(x')dx' , uj{x)=lo{V{x)) . (2.24) 
Jo 

2.3 Neutrino oscillograms of the Earth 

For given values of jAm^j, sin 2 26\3, and the type of neutrino mass hierarchy, the oscillation 
probabilities depend on the neutrino energy E v and the nadir angle of its trajectory Q u . 
Therefore a complete description of the oscillation pattern can be given by contours of equal 
oscillation probabilities in the (E u , 0„) plane. We call the resulting figures the neutrino 
oscillograms of the Earth. The plots of this type were produced for the first time by 
P. Lipari in 1998 (unpublished) and then appeared in several later publications [19, 8, 36]. 
Here we will employ this kind of plots as the main tool of our study. 

We use in our calculations the matter density distribution inside the Earth, p, as given 
by the PREM model [39] . It exhibits a characteristic structure with a relatively slow density 
change within the mantle and within the core and a sudden jump of density by about a 
factor of two at their border. Smaller jumps appear between the inner core and the outer 
one and near the surface of the Earth. The number of electrons per nucleon, Y e = N e m,N j ' p 
where m^r is the nucleon mass, equals Y e = 0.497 in the mantle and Y e = 0.468 in the 
core. For the energies E v > (1 — 2) GeV the oscillation length in vacuum l v exceeds 
1000 km, and therefore the effects of smaller-scale density perturbations are averaged out. 
Furthermore, the density profile experienced by neutrinos along any trajectory inside the 
Earth is symmetric with respect to the midpoint of the trajectory, 2 so that 

V(x) = V(L — x), (2.25) 

where L = 2i?cos @ u is the length of the trajectory. As we will see below, to a large extent 
this feature determines the properties of the oscillation probabilities and the structure of 
the oscillograms. 

2 We neglect possible short-scale inhomogeneities of the matter density distribution in the Earth which 
may violate this symmetry. 
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For numerical calculations we use the current best-fit value of Am^ = 2.5 x 10~ 3 eV; in 
the case of negligible effects of 1-2 splitting, changing Am^ is equivalent to correspondingly 
rescaling the neutrino energy. 

In Fig. 1 we present the oscillograms for the transition probability Pa for normal mass 
hierarchy and several values of sin 2 29 is. The oscillograms have two regions, separated by 
@ u = 33.1° that corresponds to the nadir angle of the border between the core and the 
mantle. For Q u < 33.1°, the neutrino trajectories cross the core of the Earth, so that 
neutrinos traverse two mantle layers and one core layer. For brevity we call this part the 
core domain of the oscillogram. Conversely, the region Q u > 33.1° corresponds to the 
mantle-only crossing trajectories, and we call it the mantle domain. As follows from the 
figure, there are several salient, generic features of the oscillation picture: 

• the MSW resonance pattern (resonance enhancement of the oscillations) for mantle- 
only crossing trajectories, with the main peak at E v ~ (5 — 7) GeV; 

• three parametric resonance ridges in the core domain, at E v > 3 GeV; 

• the MSW resonance pattern in the core domain, E v < 3 GeV, with the core resonance 
ridge at E v = 2.5 - 2.7 GeV; 

• regular oscillatory pattern for low energies: valleys of zero probability and ridges in the 
mantle domain and more complicated pattern with local maxima and saddle points in 
the core domain. 

The small windings of the contours at @ u = 65° and 70° correspond to the borders of 
the transition zone between the inner and upper mantle, while the windings at 6„ = 11° 
are due to the border between the inner and outer core. In what follows we will show that 
the whole this, apparently complicated, picture can be understood in terms of different 
realizations of just two conditions: the amplitude condition and the phase condition. 

2.4 Constant-density layers approximations 

As we shall see, the main features of the oscillograms can be well understood using the 
approximation of constant-density layers for the density distribution inside the Earth. 
Furthermore, this consideration allows us to evaluate the sensitivity of the oscillograms 
to changes of the profile. In Figs. 2 and 3 we show the oscillograms computed for two 
different approximations of this kind: 

• Approximation of fixed constant-density layers, when the Earth density profile is de- 
scribed by the core and mantle layers of constant densities that are the same for all 
neutrino trajectories. This approximation has been widely used in the literature. For 
definiteness, we have taken pi = 5.5 g/cm 3 and p2 = 11.5 g/cm 3 , which roughly cor- 
respond to the average densities in the core and in the mantle. In Fig. 2 we compare 
the PA-oscillograms calculated in this approximation and the exact results. The ap- 
proximation reproduces the oscillation pattern qualitatively well, with all the features 
present. However, quantitatively its accuracy is not high above E v > 3 GeV, in par- 
ticular, in the resonance region, and for deep trajectories B„ < 60°. For instance, one 
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0° 15° 30° 45° 60° 75° 90' 



Figure 1: Neutrino oscillograms of the Earth. Shown are the contours of constant probability Pa 
(edges of the shadowed regions) in the plane of the nadir angle of the neutrino trajectory 9„ and 
neutrino energy E v , for different values of the mixing angle 6 '13. 
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can observe an upward shift of the approximate contours compared to the exact ones 
by about 1 GeV at E u ~ (6 — 10) GeV. The shift decreases with increasing Q v . 

• Approximation of the path-dependent (trajectory-dependent) constant density layers. 
For each trajectory we find the average potentials V\ and V2 in the mantle and in the 
core, respectively, 



where Li(Q u ) is the trajectory length in the z-th layer, and then use the oscillation 
probability formulas for one or three layers of constant densities. In Fig. 3 we show the 
Pa oscillograms calculated in this approximation. One can see that the accuracy of 
the approximation is noticeably better than that of the fixed-density approximation, 
especially for the core-crossing trajectories. Again, the largest difference appears for 
the deep mantle trajectories and E v > 5 GeV. 

As can be seen from Figs. 2 and 3, the accuracy of the constant-density layers ap- 
proximations improves with decreasing neutrino energy, the reason being that the matter 
effects become less important for small E v . These approximations also in general work 
better for the core-crossing trajectories because they are dominated by the evolution in 
the core and because inside the core the density changes only by about 30%, whereas in 
the mantle it changes by about a factor of two. For mantle-only crossing trajectories, the 
approximations work better for shorter trajectories, closer to the surface. Here again the 
matter effects becomes weaker due to the vacuum mimicking phenomenon [1, 5]. 

Thus, the approximations of constant-density layers reproduce correctly the qualitative 
structure of the oscillograms obtained with the realistic PREM density profile of the Earth. 
All the resonance peaks, ridges and other structures present for the PREM-profile appear 
also in the approximate profile calculations, though their location and shape is not always 
well reproduced. Therefore one can use these approximations for understanding the main 
qualitative features of the results for realistic profiles. 

2.5 Graphical representation 

For illustrative purposes we will also use the graphical representation of the 2v oscillations 
based on their analogy with spin precession in a magnetic field [3, 40]. Let us introduce 
the neutrino "spin" vector in the flavor space s = {sx, sy, sz} with the components 



These components are essentially the elements of the density matrix. The evolution equa- 
tion for the vector s(x) can be obtained from the evolution equation for S(x), given in 




(2.26) 




Eq. (2.3): 



— = 2 B(x) x s(x) , 
ax 



(2.28) 
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Figure 2: The Pa oscillograms for the PREM density profile (colored regions; grayscale on black- 
and-white printouts) and for the fixed constant-density layers approximation (black curves). 



where 

B{x) = B(V(x)) = 4^- I sin 20, 0, 2E » V( f^ _ C os26»l . (2.29) 
4E V I Am 2 J 

In this representation the oscillations in a medium with constant density is equivalent to 
a precession of the vector s on the surface of the cone (which we will call the precession 
cone) with the axis B(V) and the opening angle 20 m . According to Eq. (2.29) the axis of 
the cone is located in the (X, Z) plane and the angle between the cone axis and the Z axis 
equals 29 m . In a medium with non-constant density profile the cone axis turns following 
the change of the angle 29 m (V(x)). The opening angle of the cone does not change in the 
adiabatic case, and it changes if the adiabaticity is broken. 

3. Physics interpretation of the oscillograms 

In this section we give a physics interpretation and description of the four main structures 
of the oscillograms mentioned in Sec. 2.3, as well as of the contours of zero probability. 

3.1 Resonance enhancement of oscillations in the mantle 

The oscillation pattern in the mantle is determined by the resonance enhancement of the 
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Figure 3: The Pa oscillograms for the PREM density profile (colored regions; grayscale on black- 
and-white printouts) and for the path-dependent constant-density layers approximation (black 
curves). 



oscillations [1, 2]. In the constant-density approximation, for a given @ u the probability Pa 
is an oscillatory function of energy, which is inscribed in the resonance curve sin 2 2Q m (E v ). 
The position of the maximum of the resonance peak is given by the MSW resonance 
condition: 

Am?, cos 20i s 

s, = E„(e,)= 2 % (e>) » , (3.i) 

where V\{Q V ) is the average value of the potential along the trajectory characterized by @j, 
(see Eq. (2.26)). Eq. (3.1) determines the resonance line in the (E u , Q u ) plane. With the 
increase of Q u the average potential decreases and consequently Er increases. According 
to Fig. 1, for the MSW resonance in the mantle we have En ~ 6 GeV. The resonance 
width AE U /Eji ~ 2tan2#i3. The condition E v = En(Q u ) ensures that the amplitude of 
oscillations is maximal and we will call it the amplitude (or resonance) condition. 

Another condition that should be met to obtain the absolute maximum of the transition 
probability, Pa = 1, is the phase condition: 

2(j)(E u , 9„) = 2u(V, E u )L(O u ) = (2k + 1)tt (3.2) 
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which means that the oscillation phase should be an odd integer of it. Since the size of 
the Earth is comparable with the refraction length in the mantle Iq = 2ir/V, the condition 
of the absolute maximum (3.2) can only be fulfilled for k = 0, i.e. when the length of 
the neutrino trajectory coincides with the half of the oscillation length in matter. The 
phase condition (3.2) gives another curve in the (E u , Q v ) plane. The intersection of the 
resonance line (3.1) and the phase condition line (3.2) gives the absolute maximum of Pa- 
Combining (3.1) and (3.2) we conclude that the absolute maximum corresponds to the 
situation when for the resonance energy the oscillation phase is tt: 

2u(V,E r )L(Q„) = tt. (3.3) 

Since L = 2RcosQ u and at the resonance 2uj = sm20is(Aml 1 /2E u ), the phase condi- 
tion (3.3) yields 

cos Q u = 9 . (3.4) 

i?sin26»i 3 Am| 1 v ; 

With the increase of #13 the peak shifts to smaller values of cosG^ (see Fig. 1). 

Let us now reformulate the resonance and the phase conditions for the case of varying 
density within the mantle. For this we first express these conditions in terms of the elements 
of the evolution matrix using the explicit results for matter of constant density, Eqs. (2.19) 
and (2.20), and then use the obtained conditions also in the case of varying matter density. 

As concerns to the resonance condition, cos 29 m = 0, it can be written according to 
Eqs. (2.19) and (2.22) as a = a*, i.e. 

&11 = $22 ) (3-5) 

or equivalently, 

ImSg^O, (3.6) 

where the superscript indicates the number of layers. This generalization, however, goes 
beyond the original MSW-resonance condition. For the constant density case it gives 

cos26» m sin0 = 0. (3.7) 

This condition has two realizations: the original one 

cos 26 m = , (3.8) 

and 

( j ) = 7 rk, (A; = 1,2,...) (3.9) 

and the latter corresponds to Pa = 0, which is realized at low energies. In the case of 
varying density (e.g., in the mantle of the Earth), there is no factorization of the reso- 
nance condition (3.6) into path-dependent and path-independent factors, as in Eq. (3.7). 
However, when the density varies slowly enough along the neutrino path, in certain limits 
Eq. (3.6) is still approximately realized in the form of Eqs. (3.8) or (3.9). At low energies 
the MSW resonance condition is not fulfilled in the Earth's mantle, and therefore the only 
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possible realization of (3.6) is the one in Eq. (3.9). In contrast to this, for large energies 
and small distances Eq. (3.9) is not satisfied, and the only way Eq. (3.6) can be imple- 
mented is through the MSW resonance (3.8). In the intermediate region neither of the two 
conditions (3.8) and (3.9) is satisfied. Thus, in the resonance region Eq. (3.6) interpolates 
between the MSW resonance condition and the condition 4> = it. 

As for the phase condition, <f> = tt/2 + wk, again we can rewrite it in terms of the 
elements of the evolution matrix, and from Eqs. (2.19) and (2.22) we get 

Re a = ReS$ = 0. (3.10) 

The absolute maximum of the transition probability occurs when both conditions (3.6) 
and (3.10) are satisfied simultaneously. In this case = and Pa = 1- This situation 
corresponds to <fi = tt/2, and from Fig. 1 we see that it is realized only for sin 2 29\3 > 0.08. 



3.2 Parametric ridges and generalized resonance condition 

For core-crossing trajectories and E u > 3 GeV the oscillatory picture is characterized by 
three ridges of enhanced oscillation probability. The ridges are the curves along which the 
probability decreases most slowly from its local or absolute maximum. The lower-energy 
ridge, which we will call the ridge 11 A" , corresponds to the energies in between those of the 
MSW resonances in the core and in the mantle, E v > (3 — 6) GeV. It was interpreted in [15, 
17, 18] as being due to the parametric enhancement of neutrino oscillations. The other two 
ridges, "B" and "C" , extend above the MSW resonance energy in the mantle. They were 
also identified as the effects of the parametric enhancement of neutrino oscillations [25]. 
Here we shall further elaborate on this interpretation. 

Recall that the parametric resonance occurs in oscillating systems with varying pa- 
rameters when the rate of the parameter change is in a special correlation with the values 
of the parameters themselves. Neutrino oscillations in matter can undergo parametric 
enhancement if the length and size of the density modulation are correlated in a certain 
way with neutrino parameters [22, 23]. This enhancement is completely different from the 
MSW effect; in particular, no level crossing is required. 

An example admitting exact analytic solution, and in fact, relevant for our discussion, 
is the "castle wall" density profile [23, 17]. This is a periodic step function with one period 
consisting of two layers of widths L\ and L2 and electron number densities N\ and N2. For 
2-flavor neutrino oscillations, the evolution matrix over one period of density modulation 
can be written as 

£( 2 ) = YI-mt-X, (3.11) 

where Y and X = (X±, A2, A3) are real parameters satisfying Y 2 + X 2 = 1 and <jj are the 
Pauli matrices in the flavor space. The oscillation probability for an arbitrary number of 
layers traversed by neutrinos can be written as a product of the amplitude 

A = X? + A 2 fxl ' (3 ' 12) 
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which does not depend on the number of the layers, and an oscillating factor, which depends 
on this number [17]. The amplitude A reaches its maximum when X3 = 0, or explicitly [17] 

X 3 = -(sic 2 cos26 1 + s 2 cicos26> 2 ) = 0, (3.13) 

where 51,2 = sin 01,2 and = cos^i^- Here 201,2 are the oscillation phases acquired in 
the layers 1 and 2, and Oip are the corresponding mixing angles in matter. Note that in 
the constant-density limit {N\ = N2 or 9\ = 62 = m ) this condition reduces to Eq. (3.7) 
with <p = (pi + (f>2. For sin0 ^ 1 it coincides with the MSW resonance condition, which is 
the maximum amplitude condition for oscillations in a matter of constant density. 

As was pointed out above, the Earth density profile seen by neutrinos with core- 
crossing trajectories can be very well approximated by three layers of constant densities, 
which is nothing but a piece of the castle wall profile; the layers "1" and "2" have to be 
identified with the Earth's mantle and core, respectively. It has been demonstrated that 
the parametric resonance condition (3.13) can be satisfied for oscillations of core-crossing 
neutrinos for a rather wide range of nadir angles both at intermediate energies [15, 16, 17] 
and high energies [25], leading to significant enhancement of the oscillation probability. 
Thus, despite the fact that the Earth's density profile consists of only three layers, the 
effects can be quite substantial. This show the relevance of the parametric resonance 
interpretation of this phenomenon, even for just "1.5 periods" of density modulation [21]. 

The parametric resonance condition (3.13) can be readily generalized to the case of 
non-constant densities in the mantle and the core of the Earth, though the generalization 
is not unique. Indeed, according to (3.11) the condition X% = can be written in terms of 
elements of the evolution matrix for two layers as the equality of the diagonal elements: 

sg> = S™. (3.14) 

We will use this equality for an arbitrary density distribution within the layers and call it 
the generalized resonance condition. Note that parameterization (3.11) of the matrix 
implies that the generalized resonance condition can also be formulated as the requirement 
that the diagonal element of be real: ImS^ = 0. 

In the right panels of Fig. 4 we show the curves of the generalized parametric resonance 
condition. Apparently, all the three parametric ridges are very well described by these 
curves. 

Since condition (3.13) does not depend on the number of layers traversed by neutrinos, 
it only ensures that the parametric oscillations occur with the maximal depth, that is, after 
some number of periods of density modulation P ~ 1. In general, the condition of maximal 
transition probability for a given number of layers need not coincide with the condition 
of maximal amplitude of the parametric oscillations. Thus, the fact that Eq. (3.13) and 
the generalized condition (3.14) give the correct description of ridges for the Earth density 
profile with three layers is rather non-trivial. 

3.3 Parametric effects and collinearity condition 

Another approach to the generalization of the parametric resonance condition to the case of 
the layers of non-constant density is based on the consideration of the evolution amplitudes 
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in the individual layers. For the particular case of two layers of constant densities, similar 
considerations have been presented in [20]. 

For density profiles consisting of two layers we have 

/ c( 2 ) c(2)\ 

S^ = S 2 S 1= J 2h b Jl) , (3.15) 




o(2)* c( 2 )* 

where 

fig) = a 2 ai - ft ft* , 5$ = °2 ft + ft a\ , (3.16) 
and ai, Pi for each layer have been defined in Eq. (2.22). The sum of the two complex 

(2) 

numbers in the transition amplitude S\ 2 can potentially lead to the largest possible result 
(if they add in the same phase and not in the anti-phase) if the two contributions to S± 2 
have the same complex phase (modulo n): 

arg(a 2 ft) = arg(/3 2 «t) mod(7r) . (3.17) 

It can also be rewritten as 

arg(aia 2 ft) = arg(ft) mod(7r) . (3.18) 

We shall call this condition the collinearity condition. It is an extremality condition for the 
two-layer transition probability under the constraint of fixed transition probabilities in the 
individual layers. In other words, if the absolute values |ft| of the transition amplitudes 
are fixed while their arguments are allowed to vary, the transition probability reaches an 
extremum when these arguments satisfy (3.17) or (3.18). For a realistic situation (neutrino 
oscillations in the Earth), changes of E v and Q u produce correlated changes of the argu- 
ments and the absolute values of the individual amplitudes, and therefore in general the 
condition (3.17) may not correspond to extrema precisely. 

If the layers 1 and 2 have constant densities, then a« = a +i cos 29i Sj, ft = — i sin 26i Si 
(see (2.22)), and the condition in Eq. (3.17) reproduces the parametric resonance condi- 
tion (3.13), i.e. X 3 = 0. Under this condition the diagonal elements of are real; 
therefore, in the case of constant-density layers the collinearity condition, the generalized 
resonance condition and the parametric resonance condition coincide. 

Denoting by Xai arm X(3i the arguments of the complex amplitudes aj and Pi, respec- 
tively, one can rewrite the collinearity condition as 

Xai + Xa2 = X/32 ~ Xf3i mod(vr) . (3.19) 

Consider now the case of three layers of in general varying densities. For the elements 
of the evolution matrix one obtains 

S[f = a 3 S$ - p 3 S[f = a 3 a 2 a x - a 3 p 2 # - ft a* 2 ft - p 3 p* 2 a x , (3.20) 
Sg } = a 3 Sf 2 ] + p 3 Sg>* = a 3 a 2 p 1 + a 3 p 2 a\ + p 3 a* 2 a\ - p 3 P* 2 ft . (3.21) 

In the case of neutrino oscillations in the Earth, the third layer is just the second mantle 
layer, and its density profile is the reverse of that of the first layer. The evolution matrix 
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for the third layer is therefore the transpose of that for the first one [38], i.e. 03 = a±, 

(3) 

/?3 = — (3 ± , and the expression for S\ 2 can be written as 

S[f = aia 2 A - ataZft + |«i| 2 /3 2 + |/?i| 2 /3 2 * . (3.22) 

Note that is pure imaginary because the core density profile is symmetric. Therefore 
the amplitude in Eq. (3.22) is also pure imaginary, as it must be because the overall 
density profile of the Earth is symmetric as well. It is easy to see that if the collinearity 
condition for two layers (3.18) is satisfied, then not only the full amplitude sf^ , but also 
each of the four terms on the right hand side of Eq. (3.22) is pure imaginary. We therefore 
conclude that if the collinearity condition is satisfied for two layers, then it is automatically 
satisfied for three layers as well. It should be stressed that this is a consequence of the facts 
that the density profile of the third layer is the reverse of that of the first layer and that the 
second layer has a symmetric profile. Once again, the collinearity of all the contributions 
to the transition amplitude potentially leads to the maximal total transition probability 
for given transition probabilities in each layer. 

Let us now confront the generalized resonance condition and the collinearity condition. 
For two layers the generalized resonance condition can be written as 

-X 3 = Im sg } = Im(a 2 ai - = (3.23) 

or, in terms of the moduli and complex phases of the amplitudes, 

\a 2 ai\ sin(xai + Xa2) = sin(x/32 - Xpi) ■ ( 3 - 24 ) 

Comparing (3.24) and (3.19) we come to the following conclusions: 

• If the generalized resonance condition (3.24) has to be satisfied independently of the 
moduli, 3 then Xai + Xa2 = Xpi ~ X/32 = ^h, and the collinearity condition (3.19) is 
satisfied too. 

• If the collinearity condition is fulfilled, then the generalized resonance condition (3.24) 
can be rewritten as 

sin( X/3 i -X/32)[|a 2 ai| ±|/3i/fc|] =0, (3.25) 
which is satisfied provided that 

X/31 = X/32 mod(vr) . (3.26) 

This condition is, in particular, fulfilled when the density profiles in both layers are 
symmetric (e.g., in the constant-density layers case), since in that case xpi = X/32 = 
±7r/2 as a consequence of T-invariance. 

Thus, we conclude that for symmetric profiles in each of the two layers (and, in partic- 
ular, for constant-density layers) the generalized resonance condition and the collinearity 

3 The same condition is required for the collinearity condition to give extrema. 
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condition coincide. For layers of non-symmetric densities, the two conditions differ. As 
can be seen in the right panels of Fig. 4, both conditions describe the parametric enhance- 
ment ridges in the oscillograms quite well. This is a consequence of the fact that the 
matter density profiles felt by neutrinos traversing the Earth can be well approximated by 
path-dependent constant density layers. 

3.4 Extrema and saddle points 

As has been discussed above, in the constant density case the positions of maxima of the 
oscillation probability are determined by two conditions: the amplitude condition and the 
phase condition. 

For matter of non-constant density, the collinearity condition or the 2-layer generalized 
resonance condition can be considered as the generalizations of the amplitude condition. 
Taking into account that the density profile of the Earth's core is symmetric (i.e. Re P2 = 0), 
one can rewrite the collinearity condition (3.18) as 

Re(aia 2 /?i) = 0. (3.27) 

Compared with Eq. (3.18), this condition has the practical advantage of not being trivially 
satisfied for = 0, which has no physically relevant implications. For trajectories crossing 
only the mantle of the Earth, one has to set «2 = 1 hi Eq. (3.27) and the collinearity 
condition reduces to Re(ai/3i) = 0. Recall that in the limit of constant matter density 
in the mantle this latter condition reduces to the condition of maximal mixing in matter 
(barring <f>i = irk). 

The phase condition cf> = vr/2 + 7rfc of the constant-density case can be generalized 
for varying density by expressing it in terms of the elements of the evolution matrix. 
According to (2.22), for one layer and the phase <j> = ir/2 + nk the product of amplitudes 
a(3* = ± sin 26 cos 29 is real. Therefore, we can generalize the phase condition Im(a/3*) 
taking instead of a and (3 the elements of the evolution matrix (amplitudes) for an arbitrary 
profile: 

lm(S u S* 2 ) = 0. (3.28) 
Furthermore, since for symmetric density profiles S12 is pure imaginary, Eq. (3.28) gives 

ReSn = 0, (3.29) 

i.e., the phase condition is fulfilled when Su is pure imaginary. 

In the graphical representation of neutrino oscillations based on their analogy with 
spin precession in a magnetic field (see Sec. 3.6), Eq. (3.28) corresponds to the condition 
that the neutrino "spin" vector is in the (sx, sz) plane. As can be easily seen, this is 
equivalent to the requirement that the transition probability be stationary with respect to 
small variations of the total distance L traveled by neutrinos: dPA/dL = 0. 

In the right panels of Fig. 4 we plot the collinearity (amplitude) condition, the gener- 
alized resonance condition and the phase condition for two different values of sin 2 2#i3 . 

As follows from the figure, the simultaneous fulfillment of the phase and amplitude 
conditions leads not only to absolute maxima of the transition probability (P4 = 1), as in 
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Figure 4: Pa oscillograms and curves of phase and amplitude conditions for the PREM profile. 
Left panels: the MSW resonance (dashed green), Reai =0 (dashed black), Im/3i = (dashed 
red) and Re 0:2 = (dashed white). Right panels: collinearity condition (solid white), generalized 
resonance condition (solid green) and phase condition (solid black). 
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the constant density case, but also to local maxima and saddle points. This is an effect of 
the multi-layer medium. To figure out why this happens, we will use the explicit results 
for three layers of constant densities. In this case [17] 

S[f = Z-iW 3 , (3.30) 

where 

Z = 2 Cl Y-c 2 , W 3 = -(2siY cos 26>i + s 2 cos 26 2 ) , (3.31) 

and 

Y = dc 2 - sis 2 cos 2(0i - 2 ) • (3-32) 
Both Z and W3 are real and therefore the phase condition (3.29) gives Z = 0, i.e.: 

2c\Y — C2 = (phase condition). (3.33) 

Also in this case 

P A = 1 - Wl (3.34) 

The collinearity condition (3.27) reduces to sin 29iSiX 3 = for constant density layers, 
and since sin 20\ 7^ 0, we have 

S1X3 = (amplitude condition). (3.35) 

Let us analyze possible realizations of these two conditions. According to (3.33), there 
are two ways to satisfy the phase condition: (i) c\ 7^ 0, Y = c 2 /2c\ and (ii) c\ = 0, 
C2 = 0. The amplitude condition (3.35) also has two realizations: (i) A3 = 0, (s\ 7^ 0), 
and (ii) s\ = 0, (A3 / 0). As we will see below, different consistent combinations of these 
realizations lead to absolute maxima, local maxima or saddle points of Pa- 

In the non-constant density case, we will use Eq. (3.27) as the general amplitude 
condition and Eq. (3.29) as the general phase condition. Using Eq. (2.22), we generalize 
the other equalities discussed above as 

ci = -► Reai = 0, (3.36) 
c 2 = Rea 2 = 0, (3.37) 
8i = -► Im/?i = 0. (3.38) 

In the left panels of Fig. 4 we show the curves that correspond to different conditions. From 
the constant-density limit it is clear why the curves Im (3i = always coincide with some 
of the amplitude condition curves: this follows from the fact that Im/3i = is a particular 
solution of the amplitude condition. 

Let us now consider all consistent realizations of the amplitude and phase conditions, 
using the terminology of the constant density approximation: 

• A3 = (si 7^ 0) (amplitude); Y = c 2 /2ci (ci 7^ 0) (phase). Plugging this expression for 
the phase condition in (3.31), we obtain 

W 3 = — . (3.39) 
ci 
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Since X% = by the amplitude condition, we have W3 = and, consequently, from 
Eq. (3.34) Pa = 1. Thus, in the constant density layers approximation a simultaneous 
fulfillment of the phase and collinearity conditions should lead to Pa = 1, provided that 
ci / or ±1. This possibility is realized only on the ridge A, at a point where the two 
curves that correspond to the generalized amplitude and phase conditions cross. The other 
curves, depicting the conditions (3.36), (3.37) and (3.38), cannot pass through this point. 
Notice that at this crossing point the oscillation half-phases in the core and mantle differ 
from it/ 2 and it, and there is a nontrivial interplay between the phases and mixing angles 
in the parametric resonance condition (3.13). 

• Ci = 0, c 2 = (the latter equality follows from Eq. (3.33)) (the phase condition). In 
this case the amplitude condition, X3 = 0, is satisfied automatically. Using the explicit 
formulas for Y (3.32), W 3 (3.31) and Eq. (3.34), we obtain Y = ±cos2(#i - 20 2 ) and 

P A = sin 2 (4#i - 26 2 ). (3.40) 

This case corresponds to the core and mantle half-phases equal to ir/2 + irk. It has a 
simple graphical representation, when it is enough to consider neutrino "spin" vector in the 
(X, Z) plane (see Sec. 3.6 below). This representation immediately leads to the expression 
in Eq. (3.40) for the transition probability [15], and shows that this probability has a 
maximum for neutrino energies between those of the core and mantle MSW resonances 
(where cos2#i > and cos2# 2 < 0), provided that (0 2 - #1) > tt/4, and above the MSW 
resonances, where cos 2#i )2 < 0. Below the resonances and between the resonances for 
(62 — 0i) < 7r/4, Eq. (3.40) corresponds to a saddle point of the transition probability. 
This agrees with the findings in [17]. 

As follows from Fig. 4, for E v < 2.5 GeV (below the MSW resonance energies) the 
intersections of the curves Re a\ = and Re a 2 = (the analogues of c\ = c 2 = in the case 
of non-constant densities within the layers) mark the positions of the saddle points. Also 
in these points the curves that correspond to the general amplitude and phase conditions 
intersect (since both conditions are fulfilled). At high energies the condition c± = is not 
satisfied: the phase in the mantle is below it. Therefore, the maxima of the transition 
probability are not achieved. 

• s\ = (ci = ±1) (the amplitude condition), Y = 2c 2 /ci (the phase condition). The latter 
equality can be written as Y = ±2c 2 . On the other hand, from the explicit expression for 
Y (Eq. (3.32)) and for s\ = one has Y = ±c 2 . Obviously, the two expressions for Y 
are consistent only if c 2 = 0. For s\ = c 2 = we find from (3.31) W3 = ±cos2# 2 , and 
consequently, 

P A = sin 2 26> 2 . (3.41) 

This realization corresponds to the oscillation half-phase in the mantle equal to it and 
therefore to the absence of the oscillation effect in the mantle (in the approximation of 
constant-density matter). The whole oscillation effect is then due to the evolution in the 
core. The oscillation half-phase in the core is a semi-integer of tt (c 2 = 0). Eq. (3.41) 
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thus simply corresponds to the maximum oscillation probability for a given mixing angle 
in matter 02- At the intersection points of the curves depicting the conditions s± = 
and C2 = the probability takes values that correspond to the resonance enhancement of 
the oscillations in the core. For non-constant density, according to Fig. 4 (left panels), 
the intersections of the curves Im/?i = and Re 02 = (which are analogues of s\ = 
and C2 = 0) correspond to local maxima. These points lie at energies below 2.5 GeV. For 
higher energies, due to the large oscillation lengths, the oscillation half-phase in the mantle 
is smaller than ir and the condition s\ = is not fulfilled. 

Notice that along the lines Im f3\ = (or <fii ps irk) the oscillation effects correspond to the 
resonance enhancement in the core, whereas the saddle points are situated along the lines 
Reai = (or <pi k, it/2 + irk). 

3.5 Absolute minima and maxima of the transition probability 

As follows from Fig. 4, the absolute minima Pa = never appear as isolated points in 
the oscillograms, but always form continuous lines (valleys of zero probability). Such a 
property (degeneracy of minima) is lifted for non-zero Am^ [37]. This is unlike for the 
absolute maxima, such as the MSW mantle peak or the parametric resonance peak in the 
core region, where instead the value Pa = 1 is reached only at a few isolated points. This 
feature is a consequence of the symmetry of the matter density profile of the Earth, and 
can be readily understood in the following way. 

The condition for the absolute minimum, Pa = 0, or Su = 0, can be written as 

ReS 12 (L,0) = 0, ImSi2(L,0) = 0, (3.42) 

and for a generic profile the absolute minima are found as the points where the curves 
corresponding to the two conditions in (3.42) intersect. However, due to the symmetry of 
the Earth's matter density profile, the condition ReS±2(L) = is satisfied automatically 
for all values of E v and @ u . Therefore the zeros of Pa simply coincide with the contour 
curves ImSuiL) = 0. 

The absolute maxima, Pa = \Sn\ 2 = 1, are realized when |Sii| 2 = 0, or 

ReSn(L,0) = 0, ImSn(L,0) =0. (3.43) 

Since in general ReSn and ImSn are independent and non-trivial functions of E u and ©1/, 
the contours that correspond to equalities in (3.43) do not coincide. The absolute maxima 
occur only at the intersections of the these contours, which explains why such maxima are 
isolated points. 

Another interesting feature of the transition probability Pa is that its dependence on 
the distance x along a given trajectory exhibits peculiar symmetry properties for trajecto- 
ries, corresponding to the absolute minima and maxima of Pa(L). Specifically, 

• For trajectories, corresponding to the absolute minima (Pa{L) = 0), 

P A (L - x) = P A (x) . (3.44) 
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That is, Pa is symmetric with respect to the midpoint of the trajectory x = L/2: 
P A {L/2 + z) = P A {L/2 -z),z< L/2. 

• For trajectories, corresponding to the absolute maxima (Pa(L) = 1), 

Pa(L-x) = l-P A (x). (3.45) 

This implies that in the middle of the trajectory Pa(L/2) = 1/2 and the function P' = 
Pa~ 1/2 is antisymmetric with respect to the midpoint: P'(L/2 + z) = —P'(L/2 — z). 

The proof is straightforward. Due to the symmetry of the density profile with respect to 
the midpoint of the neutrino trajectory, we have for any point x on the trajectory 

S(L-x,0) = S T {L,x), (3.46) 

which is essentially a consequence of T-invariance of 2- flavour neutrino oscillations [38]. 
Then, from the definition of the evolution matrix one finds S(L,0) = S(L, x) S(x, 0), 
which can be rewritten (using the unitarity of S) as S(L,x) = S(L, 0) S(x, 0)t. Plugging 
the latter relation into (3.46), we obtain 

S(L- x,0) = S*(x,0)S(L,0) . (3.47) 

For the absolute minima, Pa = 0, the evolution matrix S(L,0) should be diagonal, 
and therefore we obtain from Eq. (3.47) 

S(L-x,0) = S*(x,0)^ L e _l^. (3.48) 

Consequently, 

P A (L -x) = \S1 2 (x,0)e-^ L \ 2 = P A (x) . (3.49) 

For the absolute maxima, Pa{L) = 1, the matrix S(L,0) should be off-diagonal (see 
Eq. (3.43)) with pure imaginary elements due to the symmetry of the matter density profile. 
Therefore, using (3.47), one finds 

S(L-x,0) = S*(x,0)^.~^j , (3.50) 

and consequently, 

P A (L - x) = \S* u (x,0)(-i)\ 2 = 1 - P A (x) . (3.51) 
Notice that, according to Fig. 4, there are no local minima with Pa / 0. 

3.6 Interpretation of the oscillation pattern for core crossing trajectories 

The analysis presented in the previous subsections allows one to give a complete physics 
interpretation of neutrino oscillations in the Earth. Essentially the whole oscillatory pat- 
tern that includes the ridges, absolute and local maxima, saddle points and zeros, can 
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Figure 5: Dependence of the probability Pa (solid blue), half- phase <j> (dashed red) and the Earth 
density (gray shade) on the position along the neutrino trajectory, for sin 2 2#i3 = 0.05 and different 
values of Q v and E v . The corresponding points of the oscillogram, together with the curves of the 
amplitude (collinearity) condition, are shown in the upper-left panel. 
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be understood on the basis of different realizations of the amplitude and phase condi- 
tions. We illustrate neutrino oscillations inside the Earth for core crossing trajectories by 
Figs. 5 and 6, and the corresponding graphical representation of the oscillations is given in 
Figs. 7b~7f. For comparison, in Fig. 7a we also show the graphical representation of the 
oscillations for a trajectory crossing only the mantle of the Earth. 

Core ridge. The core resonance ridge is located at E v ~ (2.5 — 2.8) GeV. It is of the MSW 
resonance nature, but is situated below the MSW resonance line in the core: E v < Er{Q v ), 
the reason being that the values of the oscillation half-phase are different from n/2. The 
ridge does not coincides with any curve corresponding to the phase or amplitude condition. 
At one point on the ridge there is an intersection of the curves Im/?i = (the mantle 
half-phase <fi± = ir) and the collinearity condition, as well as of the Re«2 = curve, 
corresponding to the core half-phase 4>2 = vr/2 (see Fig. 4). This crossing point corresponds 
to the local maximum with zero mantle effect and maximal oscillation amplitude in the 
core. 

As follows from Figs. 5 and 6, the main contribution to the oscillation probability 
comes from the MSW-enhanced oscillations in the core, though the mantle contribution is 
not negligible. The detailed picture depends on the value of sin 2 20\^. 

For the inner core trajectories with @ u ~ 0, the core phase 4>2 > tt/2, and it decreases 
along the ridge with Q v increasing. For small mixing angles, e.g. sin 2 26>i3 < 0.05 (see 
Fig. 5, lower panels), the contributions of the two mantle layers to Pa practically cancel 
each other, so that the probability is determined by the oscillations in the core. The 
oscillation phase in the core, 02, weakly decreases with the increase of Q u , in spite of a 
substantial decrease of the trajectory length. The phase acquired in the core can be written 
as 

</>2 = — li^flL ~ const , (3.52) 

where V(Q V ) is the value of the potential averaged along a given trajectory and L2{Q V ) = 
2R C cos Q u is the neutrino path length inside the core. The ridge is below the resonance 
and therefore the oscillation length in matter decreases sharply (for small mixing) with the 
potential and also with energy. With the increase of Q u the length of the core-trajectory L2 
decreases, but simultaneously, the oscillation length decreases since the average potential 
decreases according to the PREM profile. These two dependencies compensate each other 
in Eq. (3.52) leading to weak change of 4>2- 

At Q u = 27° (point Z3) the phase reaches (j)2 = vr/2 and furthermore = it. This cor- 
responds to pure core effect. Notice that with increase of Q v the average density decreases 
and the depth of oscillations becomes smaller. In Fig. 7b we show the graphical represen- 
tation of evolution with parameters from core-ridge (point Z3). We show the precession 
cones in the mantle and in the core in the points close to border between the mantle and 
core. Shift of the evolution trajectory from the cone surfaces is due to density change and 
violation of adiabaticity. 

For sin 2 29 13 = 0.05 another intersection of lines of the amplitude condition and the 
phase conditions Re 02 = (core half-phase tt/2) occurs at E v = 3.3 GeV and <d u = 10°. 
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Evolution for this configuration is shown in Fig. 5 point A±. Notice that this point is in 
the ridge A. Here transitions in both mantle layers are non-zero but have opposite sign 
and cancel each other. So, the whole effect due to MSW resonance enhancement in core. 

For large f-3 mixing, e.g. sin 2 26*13 = 0.125, there is substantial interplay of the core and 
mantle oscillation effects (see Fig. 6). For v ~ 0, we find <fii ~ <f>2 < tt; contributions from 
two mantle layers interfere constructively and the total mantle contribution is comparable 
with the core contribution (the later is resonantly enhanced). With increase of 0„, the 
phase 4>\ increases reaching ir, whereas 4>2 decreases. Mantle contribution decreases, and 
moreover, effects of two mantle layers start to compensate each other. At Q v = 29° the 
condition of local maximum are satisfied. 

Notice that for large sin 2 2#i3 the decrease of l m with V is weaker and therefore it can 
not compensate decrease of the trajectory length. The phase stays constant if also the 
energy decreases. Consequently, for sin 2 26\^ > 0.05 the energy of ridge line decreases with 
increase of @ u . Even stronger dependence of the energy on Q u can be seen in the case of 
constant density in the core. 

Recall that in this domain of low energies the oscillations due to 1-2 splitting become 
important [37]. 

The parametric ridges differ by the oscillation phase acquired in the core, 4>2- 

Ridge A. The phase in the core <p2 ^ ^- This ridge lies in between the core resonance 
(at Q u ~ 0°) and the mantle resonance regions. At the border between the core and the 
mantle domains of the oscillogram, Q u = 33.1°, this ridge merges with the MSW resonance 
peak in the mantle. 

This ridge is large (has the largest area of probability close to 1) for small 1-3 mixing: 
sin 2 2#i3 < 0.05. As we mentioned above, for 0„ < 10° and E u ~ 3 GeV the phase in 
the core </>2 < vr/2, the mantle effect is small and two mantle layer contributions cancel 
each other. So, here we deal with resonance enhancement of oscillations in core. With 
increase of E u and @ u along the ridge both phases 4>i and 4>2 decrease. The core and two 
mantle contributions add constructively leading to large probability. For example, for the 
point A3 (Fig 5) 7Ti « 7T2 ~ 7r/3 and maximal probability build up by three comparable 
contributions. The graphical representation of this evolution is shown Fig. 7c. For the 
core part we show in some panels two precession cones, corresponding to the beginning 
and the end of the core section. The difference of these two cones reflects the effects of the 
non-constant density distribution. 

For large 1-3 mixing sin 2 2#i3 > 0.075, the region of maximal transition shifts to the 
mantle domain. Furthermore, the saddle point appear between the MSW resonance and 
the parametric resonance regions. In the saddle point A\ (Fig. 6) the phases in the core 
and the mantle 4>i ~ 02 ~ vr/2, correspond to maximal oscillatory factors. However the 
contributions from the two mantle layers have opposite sign and cancel each other. The 
corresponding graphical representation is shown in Fig. 7f. 

With increase of Q u the phases 4>i and (f)2 decrease and cancellation becomes weaker. 
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sin 2 26 = 0.125 1 




Figure 7: Graphical representation of neutrino evolution in the Earth. Shown are the trajectories 
of the neutrino vector s in the (X, Y, Z) space (solid lines) , for different values of #13 and different 
points in the oscillograms. The different panels correspond to a) the MSW peak in mantle; b) the 
point Z2 in Fig. 6; c) the point A3 in Fig. 5; d) the point B2 in Fig. 6; e) the point C\ in Fig. 6; 
f) near saddle point A\ in Fig. 6. We also show the precession cones (dotted lines) in the mantle 
and core points close to the border between the mantle and core. Sections 1-2, 2-3 and 3-4 indicate 
evolution in the first mantle layer, the core and the second mantle layer correspondingly. 
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In the ridge A for small 1-3 mixing the lines of the phase and amplitude conditions 
almost coincide that ensures stability of enhancement in large area. This stability can 
be understood also in the following way. With increasing Q„, the core segment of the 
trajectory becomes shorter and the mantle ones become longer. This is compensated by an 
increase of energy along the ridge. Indeed, since the ridge A is between the core and the 
mantle MSW resonance energies, with an increase of energy the oscillation length in the 
core (above the MSW resonance) decreases, whereas the oscillation length in the mantle 
(below the mantle MSW resonance) increases, so that the oscillation phases change rather 
weakly. Along the ridge, in the direction of larger 6„ the depth of oscillations in the core 
decreases, and this is compensated by an increase of the depth in the mantle. All this 
produces a long ridge and ensures its stability. 

Ridge B. This ridge is situated at E v > 5 GeV. For the smallest energies in the ridge 
and Q u ~ the half-phase in the core (p2 ~ (1-2 — 1.3)7r, so that oscillations in the core 
give substantial contribution effect (see Figs. 5 and 6). In the range E v = (5 — 6) GeV the 
parametric enhancement of oscillations with significant interplay of the core and mantle 
oscillation effects is realized. In Fig. 7d we show graphical representation for the point E>2 
(Fig. 6). 

With increase of Q v and E v the phase 4>2 decreases due to decrease of the trajectory 
length L2 (notice that here l m decreases too approaching Iq but much weaker than L2) In 
contrast, the phase in the mantle increases. For E v = (6 — 7) GeV one finds <p2 ~ ^ an the 
core effect is absent, so that transition occurs due to the oscillations in mantle only. 

For E v > 20 GeV the equality $2 ~ 4>i ~ vr/2 is realized. Indeed, the lines Re«i and 
Fleets as well as lines of the amplitude and phase conditions cross in the point E v ~ 21 
GeV and Q u = 31°. Here maximal enhancement of the transition probability for given 
values of mixing angles occurs [25]. It is described by Eq. (3.40). 

With an increase of both E v and Q u along the ridge the oscillation amplitudes decrease 
both in the core and in the mantle layers. The core-oscillation depth is suppressed stronger 
and therefore the region of relatively high probability lies in the same energy domain as 
the region of high probability in the mantle. 

In this ridge a qualitative picture of oscillations does not depend on sin 2 26*13. 

Ridge C. The ridge is located at E v > 11 GeV in the matter dominated region where 
mixing and consequently oscillation depth are suppressed. For Q v ~ the half-phase in the 
core equals (j>2 ~ 1.8-7T. Here main contribution to probability is due to oscillations in mantle 
whereas the core gives small (negative) contribution. We show graphical representation of 
evolution that corresponds to the point C\ in Fig. 7e. 

With increase of Q u and E v along the ridge the phase ^2 decreases and the E v > 100 
GeV reaches 37r/2, whereas the mantle phase is about (f>\ ~ ir/2. Here, as in the case of 
ridge B, the parametric enhancement of oscillations is realized [25] as described by the 
Eq. (3.40). 

Notice that qualitative features of the high energy part of ridge B and ridge C do not 
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depend on #13 since both are in the #13 factorization region where Pa oc sin 2 29\z (see 
Sec. 5.1). 



4. Approximate analytic description of neutrino oscillations in the Earth 

In this section we develop an approximate analytic approach to 2-flavor neutrino oscillations 
in the Earth. To this end, we employ a perturbation theory in the deviation of the density 
profile from that represented by layers of constant densities. This approach has been first 
developed in Ref. [30] for the description of the solar neutrino oscillations in the earth. 
Here we apply it to the study of atmospheric neutrinos. The approximation turns out to 
work extremely well in spite of the fact that the variations of the neutrino potential inside 
the Earth layers can be large | AV|/V" ~ 0.3. The high accuracy of our approach is related 
to the symmetry of the density profile of the Earth. 

4.1 Perturbation theory in AV 

Let us first consider the case of one layer of relatively weakly varying density and represent 
the matter-induced potential of neutrinos V(x) along a given trajectory as the sum of a 
constant term V and a perturbation AV(x): 

V(x) = V + AV(x) . (4.1) 

Correspondingly, the Hamiltonian of the system can be written as the sum of two terms: 

H(x) = H + AH(x) , (4.2) 

where 

/-co S 2« 8 i»2«\ AHs ^m(i »\ 

^ sin29 cos29J ' 2 \0 -l) y ' 

Here 9 = 9 m (V) is the mixing angle in matter and Q = to(V) is half of the energy splitting 
(half-frequency) in matter, both with the average potential V. Throughout this chapter 
we will denote by S(x) the evolution matrix of the system for the constant density case 
H(x) = H. The explicit expression for S(x) is given by Eq. (2.20) with 9 m = 9. 

For matter of varying density, we seek the solution of the evolution equation (2.3) in 
the form 

S(x) = S(x) + AS(x), AS(x) = -iS(x)K 1 (x), (4.4) 

where K±(x) satisfies |ifi(x) a fc| <C 1. Inserting Eq. (4.4) into Eq. (2.3), we find the following 
equation for K\(x) to the first order in AH(x) and K\(x): 

^ = St W AH W S(x) = ^ { - cos 21 (~ C ° S 2 f Si » 2 f) 
dx \ J \ J \ J 2 \ V sin2 ^ cos2 V 

+ sin 29 ( Sin2 ^- C ° S2 ^ I cos20(x) + sin 20 | ° ~M S m2<P(x)\ . (4.5) 
I cos 29 -sin20 ) vy ' i ^ V ' \ K 1 
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The first term does not contribute to S = S(L) since (Ay) = J AV(x)dx = 0, and Eq. (4.5) 
can be immediately integrated: 



Ki(L) = -sm2e\( S[n2 i C ° s2 f ] [ L AV(x) cos24>{x)dx 
y ' 2 I I cos 29 -sin2(9 J y ' YK ' 



-i 



+ . / AV(x) sm2<p(x)dx} . (4.6) 



i 



'o 



It is convenient to introduce the new variable z = x — L/2, which measures the distance 
from the midpoint of the neutrino trajectory. Then from (4.6) we obtain 

. ~ . . „ ^ I / sin 29 cos 29\ . r ( — i\ . T I , A „. 

AS = AS(L) = -t sm29 { [ - - A/+ A J > , 4.7) 

v ' I I cos2^ -sin20 ) j | ' K ' 



where 

1 fL/2 A fL/2 
AI=- AV(z) cos(2lDz) dz , A J = - / AF(z) sin(2LJz) . (4.8) 

2 J-L/2 2 J-L/2 



In these integrals, AF(z) = Al/(x(z)) and x(z) = z — L/2. Obviously, AJ vanishes 
if the perturbation AV{z) is symmetric with respect to the midpoint of the trajectory. 
Analogously, AI vanishes if AV(z) is antisymmetric. The expression for S defined in 
Eq. (4.4) with AS* given in Eqs. (4.7) and (4.8) is equivalent to Eqs. (13-16) obtained in 
Ref. [30] in the context of solar neutrino oscillations. 

Let us now consider the issue of the unitarity of the obtained evolution matrix that 
can be important for numerical calculations. Since S = S(L) is unitary, we have 

S^S = I + AS^S + 5 f AS + AS f AS. (4.9) 

The second and third terms on the RHS of this equality cancel each other thus ensuring the 
unitarity of S at the first order in AS". In order to prove this it is convenient to parametrize 
the evolution matrix S and the perturbation AS (4.7) as in Eq. (3.11): 

Y = cos(j), X = sin^(sin20, 0, -cos 20) , (4.10) 

AY = 0, AX = e (cos 29 cos£, sin£, sin 29 cos£) , (4.11) 

where 4> = <f>(L) and we have introduced 



e = sin29 \J AP + A J 2 , £ = arg(AI + % A J) . (4.12) 

It is easy to see that AS^S + S^AS = Y AY + X • AX = 0. The last term in Eq. (4.9), 
AS^ AS, violates the unitarity condition, however it is of order e 2 and therefore it does 
not break consistency of our approximation. However, for practical purposes it would be 
useful to have an expression for S which is exactly unitary regardless of the size of the 
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perturbation, even if for large perturbations the formalism is no longer accurate. In order 
to do this, we first rewrite Eq. (4.7) as follows: 



m m . f/sin 29 cos2(9\ , fo -i\ . A , Jin . 
AS = eS', S' = -i\[ - . - cos£+ . sm , 4.13 

I \ cos 26* - sm29 / \i 1 I 

and then we perform the following replacement in the expression for S: 

S = S + eS' — ► S = coseS + smeS' . (4.14) 

Note that both S' and S are unitary matrices, and that due to their specific form the 
combination on the right-hand-side of Eq. (4.14) is exactly unitary. In our computations 
we use the exactly unitary matrix (4.14). 

4.2 Application: mantle-only crossing trajectories 

We shall now use the formalism of the previous subsection to derive an explicit formula for 
the transition probability for mantle-only crossing neutrino trajectories. Since the density 
profile is symmetric with respect to the midpoint of the trajectory, the term AJ is absent. 
From Eqs. (2.20), (4.7) and (4.14) we immediately get 

Pa = [cose sin 29 sin + sine cos 29] 2 « sin 2 29 [sin <fi + AI cos 29] 2 , (4.15) 

where e = sin 29 AI and <p = 4>(L) = ujL. Here the first term in the square brackets 
describes oscillations in constant density matter with average potential V±. In order to 
obtain an explicit formula for AI, we approximate the matter density profile along the 
neutrino trajectory by a parabola: 

(4.16) 

The average value V\ and the coefficient V" depend only on the nadir angle ® v , and are 
shown in Fig. 8. Inserting the expression for AV(z) into Eq. (4.8) and integrating by parts, 
we obtain 

A/ = YlL m , m S . (4.17) 

The function f((f)) has the following features: for (f> — > (outer trajectories), one has 
f((fi) — > — 2 /15; for (j) = 7r/2 which corresponds to the maximum transition probability for 
a given mixing angle in matter, f((ft) ~ —0.13; the function |/(</>)| reaches its maximum, 
/((/>) « —0.31, at 4> — 1.07-7T. The function /(</>) changes its sign at 4> ~ 1.8-7T, and for large 
(j) it behaves as f((f>) ~ sin0/(^>. 

From Eqs. (4.15) and (4.17) it follows that: 

• in the zeroth approximation, the transition probability is given by the standard oscilla- 
tion formula for matter of constant density, with the oscillation amplitude determined 
by the mixing angle 9 = 9 m (Vi) and the phase 4> = uj(V\)L; 
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• the lowest-order correction to Pa vanishes at the MSW resonance, i.e. along the curve 
Er{Q u )- The correction changes its sign at the resonance, being positive below it 
[E v < Er(@ u )) and negative above it; 

• the largest corrections correspond to the trajectories with Q u ~ 37° (i.e. passing close 
to the core) and the phase ~ 7r. Indeed, /(0), IV/' | and L are all maximal for the 
deepest mantle trajectories 

The transition probability defined by Eqs. (4.15) and (4.17) can be cast into the form 
resembling formally the standard oscillation probability in vacuum or in matter of constant 
density: 

P A = sin 2 26' sin 2 0' , (4.18) 



where 
with 



sin 2 26' = (a 2 + b 2 ) sin 2 26, 0' = + arctan(a/6) (4.19) 

V"L - d> 2 — 3 

a = -i-r cos 26 , 6=1 + Vi'L - cos 26 . (4.20) 
40 2 12(j) 3 

This representation has the advantage that it factorizes the transition probability into a 

smooth function of neutrino energy and nadir angle, sin 2 26', and the oscillating factor 

sin 2 0' lying between and 1. In particular, zero probability curves in the oscillograms 

are given by the condition (j)'(Q u ,E u ) = irk. Factorization shows immediately continuous 

character of lines of zero probability. Note that the effective oscillation amplitude sin 2 26' 

is not actually a sine squared of any angle and should be simply understood as a notation 

defined in Eq. (4.19); in particular, it is not bounded from above by unity. 

4.3 Application: core-crossing trajectories 

For trajectories crossing the Earth's core, the evolution matrix can be factorized as 

S = £f[S2S u (4.21) 

where the subscripts '1' and '2' refer to the mantle (one layer) and core. Si and 52 
can be calculated using the formalism described in the previous sections. Since the core 
density profile is symmetric, the corresponding integral AJ2 vanishes, whereas AI2 can 
be calculated in full analogy with the derivation of AI in Sec. 4.2. In particular, we can 
approximate the core density profile by a parabola and obtain for AI2 an expression which 
is completely analogous to that in Eq. (4.17): AI2 = A/ (02, V^'L^). Here 02 = O2L2 is the 
phase acquired in the core layer. The expression for S2 is therefore 

f, ■ ■ I sin 2^2 cos 26?2 \ . „^ A T „„ N 

S 2 = cos e 2 S 2 - 1 sin £ 2 J • o5 ' 62 = sm 202 A/2 • 422 

V cos 26*2 — sin26>2 I 

For the neutrino trajectories that cross the Earth's core, it is convenient to approximate 
the density profile within each mantle layer by a linear function of the coordinate: 

V 1 (z) = V l + AV 1 (z), AV 1 (z) « V{~j- , (4.23) 

J-11 
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where L\ is the length of one mantle layer. The advantage of this parametrization is that 
AV%{z) is an antisymmetric function of z, and therefore AI X vanishes. For AJ X we obtain 
from Eq. (4.8): 

A, J i = V{L X (4.24) 
where 4>i = ui\L\ is the phase acquired in each mantle layer. Then S\ is given by 

e x = SLU29X AJ X . (4.25) 




cos 4>\+i cos 20\ sin <px sin 29 x (— i sin (pi — AJi) 
sin 2^i (— i sin 0i + A J x ) cos <j)\ — i cos 2#i sin 0i 



(4.26) 



In Fig. 8 we show the dependence of (V\, V/) and (V2, V^') on the nadir angle Q u . 
With these functions, one can find from Eqs. (4.17) and Eqs. (4.24) the quantities AI2 and 
AJ\ and then from Eqs. (4.22) and (4.25) the evolution matrices Si and S^- Substituting 
the results into Eq. (4.21), one obtains the evolution matrix for the whole trajectory. 

In the limit |AJi| <C | sin^i| Eq. (4.26) can be rewritten as 

Si = D X S X Dt, Dl ^[ 6 7 " J n /2 ) , Tl = arcsin (J£l\ , (4.27) 
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Figure 9: Contour plot of the probability Pa for the PREM density profile (colored regions; 
grayscale on black-and-white printouts) and for our analytic approximation including first-order 
corrections (black curves). 



and since AJi is real, D\ is a pure phase matrix. In this limit the total evolution matrix 
Eq. (4.21) takes the form 

S = SxD x S 2 Dx Si (4.28) 

where we have omitted the two outer matrices D\ are irrelevant for the calculation of 
Pa = |Si2| 2 . The two matrices D\ in Eq. (4.28) can be combined with S 2 to construct an 
effective core matrix, which takes into account all the first-order corrections (both in the 
core and in the mantle): 

S = Si S' 2 Si , S' 2 = Dx S 2 Dx . (4.29) 

An advantage of this approach is that both Sx and S' 2 are symmetric matrices, so that all 
the results derived previously in the constant-density approximation can be improved to 
take into account the first-order corrections by simply replacing S 2 — > S' 2 - The asymmetry 
of the mantle profile is effectively resolved. 

In Fig. 9 we compare the Pa oscillograms obtained by numerical calculations for the 
PREM matter profile with those obtained using Eqs. (4.21), (4.22) and (4.25). As can 
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be seen from this figure the degree of accuracy increases drastically once the corrections 
described here are included. 

An advantage of the described approximate analytic approach is that the coefficients 
V and V" depend solely on the nadir angle of the neutrino trajectory Q u . In particular, 
they do not depend on neutrino energy and on the values of Am 2 and 9. 

5. Dependence of oscillograms on 1-3 mixing, density profile and flavor 
channel 

5.1 Dependence of oscillograms on 1-3 mixing 

As can be seen in Fig. 1, with increasing sin 2 29\% the oscillation probability increases 
everywhere in the (E u , Q u ) plane. The evolution of the oscillation pattern appears as a 
"flow" of higher probability along nearly fixed curves towards larger values of Q u : The flow 
is along the curves of the phase condition (3.28) for the mantle trajectories, and along the 
curves of the collinearity condition for the core-crossing trajectories. These lines of flow 
move only weakly with #13. 

The region of sizable oscillation probability, P4 > 1/2, appears first for sin 2 2^13 ps 
0.009, at 6„ = 0° and E v = 2.8 GeV. It is located on the parametric ridge A. This large 
probability is due to the parametric enhancement of the oscillations. With increase of 
sin 2 2#i3 the region of sizable probability expands along ridge A, and at sin 2 29\^ ~ 0.025 
reaches the inner-mantle trajectories. At sin 2 2#i3 ~ 0.04 the regions of P4 > 1/2 appear 
also (at Q u = 0°) in the core MSW resonance region and in the parametric ridge B. The 
curves of zero oscillation probability change only very slightly with #13 . The nature of the 
increase of the oscillation probability with increasing #13 is different in different regions of 
the parameter space. 

Analytic expressions for various structures in the oscillograms, derived in Sec. 2.3, allow 
one to understand this evolution. Two features determine dependence of the oscillograms 
on 1-3 mixing: 

• factorization of the #13 dependent factors in the probability, and 

• dependence of the amplitude and phase conditions on 1-3 mixing. 

Let us consider first the regions outside the resonances: E v > E R + A£^ ax , and 
E v < E R -AE™ X . Here AE™ ax = E R tan 29[f x and tan 29f£ x ~ (0.1-0.2) corresponds to 
the maximal allowed values of 1-3 mixing. In practice, for the mantle-crossing trajectories 
these are the regions with E v < 5 GeV and E v > 8 GeV. In these regions, the mixing 
parameter in matter can be approximated by 



sin 2 29. 
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does not depend on 1-3 mixing. So, for the mantle domain trajectories the lines of constant 
phase that determine the lines of flow do not depend on 1-3 mixing. As follows from (5.1) 
and (5.2), the oscillation probability for one layer (mantle) factorizes: 



Pa ~ sin 2 26*13 sin 2 4> 



1 



(5.3) 



1 - 2VE l/ /Am 2 



Therefore, the probability increases uniformly in the whole this area and the lines of zero 
and maximum transition probability do not move. For the core crossing trajectories similar 
analysis holds for E v < 2 GeV. 

Let us first consider the parametric resonance condition (3.13) that determines ap- 
proximately the lines of "flow". Beyond the MSW resonance regions, the phases and 
therefore q and Sj (i = 1,2), do not depend on the 1-3 mixing. The mixing angles in 
matter enter the condition as cosines cos 26V Outside the resonance regions Oi ~ #13 <C 1 
or Oi ~ it/ 2 and therefore cos 26* j ~ ±1 weakly depends on 1-3 mixing. So, the condition 
(and therefore the lines of flow) shifts only weakly with change of #13. 

Let us now consider the resonance regions (2 < E v < 12) GeV. Here dependence of the 
transition probability on sin 2 20 13 is non-linear. Not only the depth of the oscillations, but 
also the oscillation length in matter in each layer depends on sin 2 26*13 substantially, and 
the latter influences the interference effects. As a result, the transition probability changes 
with 6*i3 differently in different regions, and in addition the "lines of flow" shift. 

• The MSW resonance peak in the mantle is determined by the condition (3.1). With 
increasing sin 2 26*13, the peak shifts to larger values of Q u and to slightly larger ener- 
gies. This can be readily understood. Indeed, the oscillation length at the resonance 
decreases with increasing 6*13 as l m = l v j sin 26*13 and therefore the condition <pi = tt/2 
is satisfied for shorter (more external) trajectories. For these trajectories the average 
density becomes smaller, so that the resonance energy increases: Er oc 1/V. The 
width of the resonance peak at half height increases as tan 26*13 both in neutrino en- 
ergy and in Q u variables (see the discussion below Eq. (3.1)). Using the resonance 
condition (3.1), we obtain from (3.4) 



where d is the column density that corresponds to maximum of the transition prob- 
ability. Eq. (5.4) gives an immediate relation between #13 and the nadir angle of the 
neutrino trajectory on which the absolute maximum of the transition probability is 
realized. This can, in principle, be used for measuring 6*13: the method would simply 
consist in the determination of <d u of the mantle-only crossing trajectory corresponding 
to the absolute maximum of the conversion probability Pa — 1 — Pee- 

• Core ridge slightly shifts with increase of 6*13 to larger energies, especially at Q v « 0. 

• For ridge A, the region of sizable transition probability and the position of the max- 
imum also move towards smaller |cos6v|. With increasing #13 the oscillation length 
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decreases, especially in the resonance region. Therefore, the same phases in the core 
and mantle can be obtained for shorter trajectories in the core and therefore for larger 

e„. 

• Ridge B evolves weaker with #13: the energy corresponding to the maximum of the 
transition probability stays rather close to that of the mantle MSW resonance, E v ~ 
(5 — 6) GeV. At higher energies the ridge is in the factorization region. With increase 
of #13 at Q u « the ridge shifts to smaller E v . 

• The situation for ridge C is similar to ridge B. 
5.2 Dependence on the Earth's density profile 

In different parts of the (E u , 0„) plane the oscillation probabilities have different sensitivity 
to the modifications of matter density profile. The sensitivity is very weak (independently 
of the form of perturbation) in the following parts: 

• Q u > 84°: here the length of the trajectory, and therefore the oscillation phase, are 
small, so that effect of "vacuum mimicking" [1, 5] takes place. To a good approximation 
matter does not affect the oscillation probabilities, irrespectively of whether or not 
the matter-induced potential V is small compared to the kinetic energy difference 

• E v < 2 GeV: here one has V <C Am^ 1 /2E V , and so matter effects on the oscillations 
driven by 1-3 mixing and splitting are small for all values of Q u . 

The region of high sensitivity to the density profile of the Earth is bounded by E v > 
(3 — 4) GeV and 0„ < 66°. This is the region of energies close to and above the resonance 
energies (matter dominance) and of sufficiently long trajectories. The latter condition 
ensures an accumulation of the matter effects over long distances and therefore a sensitivity 
to large scale structures of the density profile. For the core-crossing trajectories the border 
of the sensitivity region is lower: E v ~ 2 GeV. Changes of the oscillograms in this region 
depend on the specific form of the modification of the matter density distribution. 

For illustration, we consider here the effects of three different modifications of the 
Earth density profile: (1) replacing PREM profile by constant-density mantle and core 
layers, (2) modification of the core/mantle density ratio, and (3) changes of the position of 
the border between the mantle and the core. 

1. In Fig. 2 we show the results of the fixed constant-density layers approximation, char- 
acterized by the constant potentials V\ and V2 in the mantle and core. This profile can be 
considered as an extreme case of flattening of the mantle and core density distributions. 
For the mantle-only crossing trajectories there are two lines in the oscillograms where the 
oscillation probabilities for the two profiles are equal: (i) Q u « 53°, which corresponds to 
the trajectory with Vprem = Vi, and (ii) the resonance energy curve E v ~ Er(® v ). In- 
deed, along the resonance curve the first-order correction to Pa due to the deviation of the 
density profile from the averaged constant one disappears (see the discussion in Sec. 4.2). 
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For Q u < 53° the oscillation pattern for the fixed constant-density profile is shifted to 
higher energies compared to that for the PREM profile. Furthermore, for E v > Er(Q v ) 
one has P con st > -Pprem > whereas for lower energies, E v < E R (Q U ), P con st < Pprem, which 
essentially reflects the shift of the resonance peak when profile is changed. The shift of 
contours in the energy scale increases as Q v decreases and reaches maximum for the deepest 
mantle trajectories. Quantitatively, in the high energy region, E v > 6 GeV, for Q u ~ 37° 
the changes of probability equal P con st/PpREM = 0.2/0.06, 0.3/0.1, 0.4/0.2, 0.5/0.3, 0.9/0.5. 
That is, the modification is characterized by factor 2-3. The shift of contours in energy 
scale is about 1 GeV at E v ~ 10 GeV. The differences between the results for the two 
profiles weakly depend on 6*13. For Q v > 53°, the PREM profile gives higher probabilities 
than those for the fixed constant-density layers approximation for E u > Er(@ u ) and lower 
probability for E v < Er(@ v ). 

For core crossing trajectories the size of changes is similar. 

2. Fig. 10 illustrates the effects of the increase (decrease) of the core density: V 2 (x) — > 
kV2(x), where k = const. The shapes of the density profiles in the mantle and core 
are taken according to the PREM profile. The total mass of the Earth is unchanged 
and therefore the density of the mantle should be reduced (increased) correspondingly: 
AVl/Vl = -[(R/R c f - 1]- 1 (V 2 /V L )AV 2 /V2 « -0.3AV 2 /V 2 . The effects of the increase 
and decrease of the core potential on the oscillogram are opposite and for definiteness we 
will consider the case of an increase of V 2 , so that for mantle-only crossing trajectories V\ 
decreases. Since the resonance energy is proportional to 1/V, a decrease of density leads 
to a shift of the resonance peak to higher energies and also to smaller Q u , to satisfy the 
phase condition <p = ir/2. The strongest effect is for Q u < 60° and above the resonance. 
Quantitatively, a ~ 6% decrease of the mantle density (corresponding to the 20% increase 
of V 2 ) leads, for 6„ < 37°, to an about 1 GeV upward shift of the contour curves. As follows 
from the figure, the increase of probabilities with respect to those for the PREM profile is 
^20%/-Pprem = 0.1/0.03, 0.2/0.1,0.4/0.2, 0.5/0.3,0.7/0.5. The curve of unchanged proba- 
bility is close to the resonance curve. Below the resonance energy the PREM profile leads 
to larger probabilities than the modified one. 

Notice that these changes are rather similar to the changes for the constant-density 
layers profile. The reason is that in both cases in the inner parts of the mantle the density 
is smaller than that given by the PREM density profile. 

The effect of density modifications on the transition probability is nearly linear in the 
variable of density change. 

For core-crossing trajectories, an increase of the core density leads to a shift of the 
oscillatory pattern to larger <d u and higher energies, and also results in an increase of 
probability. These features can be understood using the parametric resonance condition. 

3. Fig. 11 illustrates the effects of increase (decrease) of the core radius. Again, the 
requirement of keeping the total mass of the Earth fixed imposes a rescaling of the overall 
core and mantle densities. As in the previous case, this rescaling is what induces the most 
important effects on the oscillogram. 
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Figure 10: Pa oscillograms for the PREM density profile (colored regions) and for a 20% variations 
of the core/mantle density ratio. The total mass of the Earth is kept fixed. 



The sensitivity of the oscillation probabilities to the variations of the matter density 
distribution in the Earth can in principle be used for studying the Earth interior with 
neutrinos, i.e. to perform an oscillation tomography of the Earth [37]. Using the features 
described above one can work out the criteria for the selection of events that are sensitive 
to a given type of variations of the profile. For instance, in the case of flattening of the 
density distribution, one can divide the whole parameter space into 4 parts: (1) @ u < 53°, 
E v > E R ; (2) 6„ > 53°, E v < E R ; (3) Q v < 53°, E v < E R ; (4) @ u > 53°, E v > E R , and 
select a sample of the edike events sensitive to these regions. The flat profile gives larger 
number of events then PREM profile in (1) and (2) regions and smaller number of events 
in the regions (3) and (4). This feature allows one to distinguish effect of flattening from 
other possible modifications of the profile as well as effect of uncertainty in #13. 

Notice that the changes of the oscillation probability due to modifications of the 
core/mantle density ratio and of the position of the border between the mantle and the 
core can be noticeable if the magnitude of the modifications is as large as 20% (see Figs. 10 
and 11). However, to detect them experimentally, the energy and nadir angle resolutions 
should be high as well as large event statistics is necessary. Therefore, the oscillation 
tomography of the Earth will put very challenging demands to future experiments. 
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Figure 11: Pa oscillograms for the PREM density profile (colored regions) and for a 20% variations 
of the core radius. The total mass of the Earth is kept fixed. 



5.3 Probabilities for other oscillation channels 

In Figs. 12 and 13 we show the oscillograms for the other channels. As follows from 
the discussion in Sec. 2.2, in the approximation of zero 1-2 splitting (AmJx = 0) all the 
probabilities that involve v e depend on the single probability Pa- This is related to the 
fact that v e is unchanged upon going from the flavor basis to the propagation one. 

The probability P ee = 1—Pa is just complementary to Pa with all the features inverted. 
P efJ , is just Pa scaled by the factor s| 3 (see the middle panel in Fig. 12). The maximal value 
of this transition probability is therefore s^. Similarly, P eT is Pa scaled by the factor c^- 

The probabilities of transitions that do not involve z/ e have more complicated structure 
since now both the initial and the final states do not belong to the propagation basis and 
therefore some additional interference occurs. The survival probability Pan in Eq. (2.16) 
can be written as 

Pau = Pvac ~ S^Pa + 2S23 C 23 ( Re ^33 ~ COS 2<j) vac ) , (5.5) 

where P vac = 1 — sin 2 2#23 sin 2 4> vac is the usual 2-flavor vacuum oscillation probability and 
4> vac = Am'^ 1 L/4E u is the vacuum oscillation phase. P vac describes the oscillation effect 
in the absence on the 1-3 mixing. The other two terms in (5.5) describe the effects of 
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Figure 12: P ee , P e/J and -P MM oscillograms, for neutrinos (left panels) and antineutrinos (right 
panels), sin 2 29 13 = 0.05, sin 2 2# 2 3 = 1 and Am^ = 0. 
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the 1-3 mixing. We can use the approximate results of Sec. 4 to estimate Re ^33. Recall 
that in Sec. 4 we used a symmetric Hamiltonian for 2v system, which differs from the one 
introduced in Sec. 2.1 by a term proportional to the unit matrix (2.18). The term Re ^33 
essentially reflects evolution of the 1-3 system with respect to the state v 2 and so the Zv 
form of the Hamiltonian should be restored. According to Eq. (2.18) the relation between 
the corresponding evolution matrices (up to the removed decoupled state) is 




Then 

^33 = S 22 , (5.7) 

where S22 is the 22 element of the matrix (2.20) with the correction given in (4.7). 
Let us consider the mantle trajectories. Then from (5.7) we get explicitly 

ReA 33 = cos(0 + V) + 2 sin 2 013 sin V> [sin 0+ 2 cos 2 % A/] , (5.8) 

where the phase <fi and the correction integral AI are defined in (2.20) and (4.8), respec- 
tively, and the angle #13 is the 1-3 mixing angle in matter calculated at the average value 
of the matter density. Then the survival probability for the mantle-only crossing 
trajectory can be written as 

= 1 - sin 2 26 23 sin 2 (^±^j - s\ 3 P A 

+ - sin 2 26*23 sin ip (2 sin 2 £13 sin <f> + sin 2 26>i 3 Al) . (5.9) 

Notice that the first two terms here correspond to 2-flavor vacuum probability with 
the modified phase. The probability can be rewritten as 

Pw = Pvac + AP (5.10) 

with 



AP = -s% 3 Pa - 2s| 3 c| 3 [ cos 2(/> vac - cos(0 + V) 

— 2 sin 2 6^13 sin ip sin <f> — sin 2 2^13 sin ip Al] . (5.11) 

In the limit S13 — > one has 

*"jf(£-^)* 

which, together with (5.6), implies <f> + tp = 2</> mc . 

These formulas can be used for the analysis of numerical results shown in Figs. 12 
and 13. The first correction term in Eq. (5.11) is negative, so that it reduces the survival 
probability. The strongest effect of this mixing occurs in the resonance region. For high 
energies the 1-3 mixing is suppressed and P^ is again described well by the vacuum 
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oscillation formula. Notice that in some regions with zero P vac the 1-3 mixing leads to 
the positive contribution which is due to the last term in (5.5) so that the lines of of zero 
probability are interrupted. Also there are no lines of maximal probability, as zero lines 
for P A . 

All the correction terms to P vac are in general of the same order; therefore the modi- 
fication of the vacuum probability is rather complex and not immediately related to P^. 
The corrections are in general large in the regions of the parameter space where P eiI is 
large, i.e. in the resonance regions and in the places of the resonance peaks and parametric 
ridges of Pa- Instead of continuous lines of maximal probability, in the case of non-zero 
1-3 mixing one obtains saddle points and local maxima. In particular, the saddle point of 
the P^n probability appears in the mantle resonance region. 

The structure of the transition probability P^ T is similar to that of P^ . 

In the antineutrino channel, the 2-flavor probability P2 is suppressed by matter and 
the strongest transitions occur in the vacuum oscillation region. Some interference effects 
are seen for the core domain but no substantial parametric enhancement is realized. 

6. Discussion and conclusions 

1. We have worked out a detailed and comprehensive description of neutrino oscillations 
driven by non-zero 1-3 mixing inside the Earth. The description is given in terms of the 
oscillograms of the Earth: contours of constant oscillation probabilities in the (E u , 0„) 
plane. In this first publication we have neglected the 1-2 mass splitting Am^, which is a 
good approximation for high neutrino energies, E v > 1 — 2 GeV. 

2. We found that the oscillograms have a regular structure with several generic features: 
(i) the MSW peak in the mantle domain of the oscillogram; (ii) the MSW peak (ridge) in 
the core domain, (iii) three parametric ridges in the core domain; (iv) regular oscillatory 
pattern at low energies that has different features in the core and in the mantle domains 
of the oscillograms. We presented a detailed description of these structures: their position 
in the (E Ul 0„) plane, their evolution with changing 1-3 mixing and their dependence on 
the density profile of the Earth. 

The most interesting features of the oscillograms appear at relatively high energies, 
E v = (2 — 12) GeV, and Q u < 75°, i.e. in the resonance region. Notice that this region is 
not covered by the existing or forthcoming accelerator experiments, and on the other hand, 
statistics in the current atmospheric neutrino experiments is too low. All these experiments 
can only study very small effects on the "tails" of those interesting oscillation phenomena. 
Exploration of the resonance regions thus constitutes a significant experimental challenge 
for future experiments. 

3. We studied the accuracy of the calculations based on the constant-density layers ap- 
proximation to the Earth matter density profile. This approximation reproduces all the 
features of the oscillograms qualitatively well, though there are some quantitative differ- 
ences and shifts of the structures in the (E u , @ v ) plane. The strongest deviations appear 
in the domain of deep mantle trajectories and high energies {E v > 6 GeV). 
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4. We showed that a complete physics interpretation of the oscillograms can be given 
in terms of different realizations of just two conditions: (i) the amplitude (or resonance) 
condition and (ii) the phase condition. In the case of one layer of constant density these 
conditions are reduced to the MSW resonance condition and half-phase equality 4> = tt/2 + 
irk. They determine the position of the absolute maxima of the transition probability 
For three layers of constant densities the amplitude condition is reduced to the parametric 
resonance condition. This condition describes the position of the parametric ridges. The 
phase condition gives the position of the maximum along the ridge. 

We show that the parametric resonance condition formulated for two layers describes 
the extrema and the ridges in the three layers case. This is a consequence of the symmetry 
of the overall matter density profile. In the case of three layers, the amplitude and the 
phase conditions determine not only the absolute maxima of the transition probability (as 
in the one-layer case), but also its local maxima and saddle points. 

5. We generalized the amplitude and phase conditions to the case of varying densities 
in each layer. The generalization is not unique. In this connection we introduced the 
generalized resonance condition and the collinearity condition. The two conditions coincide 
in the case of constant density layers with symmetric overall profiles but differ in the 
non-constant density case. We showed that both these generalized conditions describe 
the positions of various structures of the oscillograms, in particular, of the extrema, very 
accurately. 

6. We derived approximate analytic formulas for the probabilities. For this we have 
developed a perturbation theory in deviations from the constant density AV/V. We showed 
that already the first order approximation in AV/V reproduces the oscillograms for realistic 
(PREM) density profile of the Earth with a high precision. Again, the symmetry of the 
density profile plays the key role. 

7. We studied the dependence of the oscillograms on #13. The changes of the oscillograms 
with increasing #13 have a character of flow of high probabilities towards the regions of 
larger Q u . The lines of flow shift only weakly with changing #13. We found that the 
transition probability Pa can be of the order 1 for sin 2 2^13 as small as 0.01. Therefore, 
even if the next generation of reactor and accelerator experiments fail to find non-zero 1-3 
mixing, significant oscillation effects due to this mixing may still show up in atmospheric 
neutrino data. Those are expected in the region E v ~ 3 — 5 GeV and Q u pa 0° — 26°. 

8. We studied the dependence of the oscillograms on the density profile of the Earth. We 
found the regions in the {E v , Q u ) plane where the sensitivity to various perturbations of 
the density profile is maximal, and we identified the corresponding effects. In particular, 
the dependence of the oscillograms on flattening of the density distributions inside the 
layers, on the changes of the overall densities of the core and mantle and on the position 
of the border between the mantle and the core has been quantified. This analysis can be 
used for discussions of the oscillation tomography of the Earth. 
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9. The oscillograms for different flavor channels as well as for neutrinos and antineutrinos 
have been constructed and their properties discussed. 

Various applications of the results obtained in this paper will be presented in forthcoming 
publications. 
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